Random close packing revisited: How many ways can we pack frictionless disks? 
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We create collectively jammed (CJ) packings of 50-50 bidisperse mixtures of smooth disks in 2d 
using an algorithm in which we successively compress or expand soft particles and minimize the 
total energy at each step until the particles are just at contact. We focus on small systems in 2d and 
thus are able to find nearly all of the collectively jammed states at each system size. We decompose 
the probability P((j>) for obtaining a collectively jammed state at a particular packing fraction <f> 
into two composite functions: 1) the density of CJ packing fractions p(<j>), which only depends on 
geometry and 2) the frequency distribution /3(4>), which depends on the particular algorithm used to 
create them. We find that the function p{(f>) is sharply peaked and that /3((j>) depends exponentially 
on <f). We predict that in the infinite system-size limit the behavior of P((j>) in these systems is 
controlled by the density of CJ packing fractions — not the frequency distribution. These results 
suggest that the location of the peak in P((f>) when N — > oo can be used as a protocol-independent 
definition of random close packing. 

PACS numbers: 81.05.Rm, 82.70.-y, 83.80.Fg 



I. INTRODUCTION 

Developing a statistical mechanical description of 
dense granular materials, structural and colloidal glasses, 
and other jammed systems Q composed of discrete 
macroscopic grains is a difficult, long-standing problem. 
These amorphous systems possess an enormously large 
number of possible jammed configurations; however, it is 
not known with what probabilities these configurations 
occur since these systems are not in thermal equilibrium. 
The possible jammed configurations do not occur with 
equal probability — in fact, some are extremely rare and 
others are highly probable. Moreover, the likelihood that 
a given jammed configuration occurs depends on the pro- 
tocol that was used to generate it. 

Despite difficult theoretical challenges, there have been 
a number of experimental and computational studies that 
have investigated jammed configurations in a variety of 
systems. The experiments include studies of static pack- 
ings of ball bearings 0, Q , slowly shaken granular ma- 
terials 0, |j| , sedimenting colloidal suspensions |(J , and 
compressed colloidal glasses Q. The numerical studies 
include early Monte Carlo simulations of dense liquids [g , 
collision dynamics of growing hard spheres M, serial de- 
position of granular materials under gravity jloL ITTl IT^ | , 
various geometrical algorithms |l3l Il4l Il5| , compression 
and expansion of soft particles followed by energy mini- 
mization [l^, and other relaxation methods [l7| . 

The early experimental and computational studies 
found that dense amorphous packings of smooth, hard 
particles frequently possess packing fractions near ran- 
dom close packing (f> Tcp , which is approximately 0.64 in 3<i 
monodisperse systems and 0.84 in the 2d bidisperse 
systems discussed in this work [TH Il9j|. However, more 
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FIG. 1: The probability distribution P((j>) to obtain a col- 
lectively jammed state at packing fraction <f> in 2d bidisperse 
systems with N — 18 (dotted), 32 (dashed), 64 (dot-dashed), 
and 256 (solid). 

recent studies have emphasized that the packing fraction 
attained in jammed systems can depend on the process 
used to create them. Different protocols select particular 
configurations from a distribution of jammed states with 
varying degrees of positional and orientational order |20| . 

Recent studies of hard particle systems have also shown 
that different classes of jammed states exist with differ- 
ent properties 21,]. For example, in locally jammed (LJ) 
states, each particle is unable to move provided all other 
particles are held fixed; however, groups of particles can 
still move collectively. In contrast, in collectively jammed 
(CJ) states neither single particles nor groups of particles 
are free to move (excluding 'floater' particles that do not 
have any contacts). Thus, CJ states are more 'jammed' 
than LJ states. 

In this article we focus exclusively on the properties of 
collectively jammed states. These states are created us- 
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ing an energy minimization procedure 16, 19] for systems 
composed of particles that interact via soft, finite-range, 
purely repulsive, and spherically symmetric potentials. 
Energy minimization is combined with successive com- 
pressions and decompressions of the system to find states 
that cannot be further compressed without producing an 
overlap of the particles. As explained in Sec.|nl this pro- 
cedure yields collectively jammed states of the equivalent 
hard-particle system. 

In previous studies of collectively jammed states cre- 
ated using the energy minimization method, we showed 
that the probability distribution of collectively jammed 
packing fractions narrows as the system size increases 
and becomes a ^-function located at </>o in the infinite 
system-size limit 0|. We found that 4>q was simi- 
lar to values quoted previously for random close packing 
|l8j. The narrowing of the distribution of CJ packing 
fractions as the system size increases is shown in Fig. ^ 
for 2d bidisperse systems. However, it is still not clear 
why this happens. Why it is so difficult to obtain a col- 
lectively jammed state with <f> ^ O hi the large sys- 
tem limit? One possibility is that very few collectively 
jammed states exist with <p ^ 0o- Another possibility is 
that collectively jammed states do exist over a range of 
packing fractions, but only those with packing fractions 
near 0o are very highly probable. 

Below, we will address this question and other related 
problems by studying the distributions of collectively 
jammed states in small bidisperse systems in 2d. For 
such systems we we will be able to generate nearly all 
of the collectively jammed states. Enumeration of nearly 
all CJ states will allow us to decompose the probability 
density P(4>) to obtain a collectively jammed state at a 
particular packing fraction <p into two contributions 

P{<f)=p{cj > )p{<j > ). (1) 

The factor p((f>) in the above equation represents the den- 
sity of collectively jammed states (i.e., p((f>)d(f> measures 
how many distinct collectively jammed states exist within 
in a small range of packing fractions d<j)). The factor j3 
denotes the effective frequency (i.e., the counts averaged 
over a small region of (f>) with which these states occur. 

We note that the density of states p(<f>) is determined 
solely by the topological features of configurational space; 
it is thus independent of the the protocol used to gener- 
ate these states. In contrast, the quantity f}{4>) is proto- 
col dependent, because it records the average frequency 
with which a CJ state at 4> occurs for a given protocol. 
For example, for algorithms that allow partial thermal 
equilibration during compression and expansion, the fre- 
quency distributions are shifted to larger <f> compared to 
those that do not involve such equilibration. 

The decomposition Q will allow us to determine which 
contribution, p(</>) or /3(0), controls the shape of the prob- 
ability distribution P{4>) in the large system limit. Oth- 
ers have studied the inherent structures of hard-sphere 
liquids a nd glass es, but have not addressed this specific 
question |22ll23l|. We will show below that p(<f>) controls 
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FIG. 2: Five distinct collectively jammed states that exist 
at the same packing fraction <j) = 0.81073 for a 2d bidisperse 
system with N = 12 particles. In (a) through (d), unshaded 
particles are in the same positions, while the shaded particles 
are in different locations from panel to panel. Particles labeled 
1, 2, and 3 are in different positions in panels (d) and (e), while 
the other particles are in the same positions. 

the width of the distribution of CJ states in the infinite 
system-size limit. We also have some evidence that the 
location of the peak in P((j>) in the large N limit is also 
determined by the large N behavior of p(<p) . We will also 
argue that for many procedures the protocol-dependence 
of the frequency distribution /3(<f>) is too weak to sub- 
stantially shift the peak in P(4>) for large systems. Thus, 
our results suggest that for a large class of algorithms the 
location of the peak in P(<fi) can be used as a protocol- 
independent definition of random close packing in the 
infinite system size limit. 



II. METHODS 

Our goal is to enumerate the collectively jammed con- 
figurations in 2d bidisperse systems composed of smooth, 
repulsive disks. We will focus on bidisperse mixtures 
composed of N/2 large and N/2 small particles with a 
diameter ratio r) = 1.4 because it has been shown that 
the se sy stems do not easily crystallize or phase separate 
[151 Hq . We consider system sizes in the range N = 4 to 
256 particles. For N < 10, we were able to find nearly 
all of the collectively jammed states. For N = 12 (14) we 
found more than 90% (60%) of the total number. Since 
the number of collectively jammed states grows so rapidly 
with N, we are not able to calculate a large fraction of 
the CJ states for N > 14, but as we will show below, we 
can still make strong conclusions about the shape of the 
distribution of CJ states in large systems. 

We utilize an energy-minimization procedure to cre- 
ate collectively jammed states [l^. We assume that the 
particles interact via the purely repulsive linear spring 
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FIG. 3: Fraction /a of distinct collectively jammed states 
with an excess number of contacts A = N c — iV" lln over a range 
of system sizes, N = 6 (circles), 8 (squares), 10 (diamonds), 12 
(upward triangles), and 14 (downward triangles). For N = 6, 
all CJ states have A = 0. 

potential 

V(rij) = |(1 - r y -/^-) 2 e(^/r y - - 1), (2) 

where e is the characteristic energy scale, is the sep- 
aration of particles i and j, = {di + dj) /2 is their 
average diameter, and 0(x) is the Heaviside step func- 
tion. The potential J5J, is nonzero only for r\j < dij, 
i.e., when the particles overlap. Jammed states are ob- 
tained by successively growing or shrinking particles fol- 
lowed by relaxation via potential energy minimization 
until all particles (excluding floaters) in the system are 
just at contact. In these prior studies, we showed that 
the distribution of collectively jammed states does not 
depend sensitively on the shape of the repulsive poten- 
tial V(rij). Note that our process for creating jammed 
states differs from the fixed volume energy minimization 
procedure implemented in Ref. 16]. In the description 
below, the energies and lengths are measured in units of 
e and the diameter of the smaller particle d\ . 

For each independent trial, the procedure begins by 
choosing a random configuration of N particles at an 
initial packing fraction <f>i in a square box with unit length 
and periodic boundary conditions. The positions of the 
centers of the particles are uncorrelated and distributed 
uniformly in the box. We have found that the results do 
not depend on the initial volume fraction fa as long it is 
significantly below the peak in p((f>). We chose fa = 0.60 
for most system sizes. 

After initializing the systems, we find the nearest local 
potential energy minimum using the conjugate gradient 
algorithm |24| . We terminate the energy minimization 
procedure when either of the following two conditions is 
satisfied: 1) two successive conjugate gradient steps n 
and n + 1 yield nearly the same total potential energy 
per particle, (V n+1 - V n )/V n < 5 = 10~ 16 or 2) the total 
potential energy per particle is extremely small, V n +i < 
V ■ — 10 -16 

Following the potential energy minimization, we decide 
whether the system should be compressed or expanded to 
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FIG. 4: Fraction of CJ states n a /ri}'° t versus the ratio of the 
number of trials nt to the total number of trials n l ° l required 
to find all CJ states nl ot for several system sizes: N = 6 
(dot-dashed), 8 (dashed), 10 (thick solid), 12 (long dashed), 
and 14 (dotted). The curves for N = 10, 12, and 14 collapse. 
The thin solid line is a least-squares fit to Eq. ©. The dia- 
mond, upward triangle, and downward triangle symbols give 
the maximum number of trials attempted for N = 10, 12, and 
14, respectively. 

find the jamming threshold. If V n +i > V max — 2 x 10~ 16 , 
particles have nonzero overlap and thus small and large 
particles are reduced in size by Adi — d\AfaH2fa) and 
Ag?2 = J/Adi, respectively. If, on the other hand, 
Vn+i < VJ n i n , the system is below the jamming thresh- 
old and all particles are thus increased in size. After the 
system has been expanded or compressed, it is relaxed 
using potential energy minimization and the process is 
repeated. Each time the procedure switches from expan- 
sion to contraction or vice versa, the packing fraction 
increment Acj) is reduced by a factor of 2. The initial 
expansion rate was Afa = 10 . 

When the total potential energy per particle falls 
within the range F ma3! > V > V m in> the process is termi- 
nated and the 'jammed' packing fraction is recorded. If 
the final state contains floater particles with 2 or fewer 
contacts, we remove them, minimize the total potential 
energy, and slightly compress or expand the remaining 
particles to find the jamming threshold. Note that the 
final configurations are slightly compressed with overlaps 
in the range 10 -9 < 1 — ry/dy < 10 -8 . We have verified 
that our results do not depend strongly on the parame- 
ters Kiin, Knax, and Afa. 

For each system-size N, this process is repeated us- 
ing n' t independent random initial conditions and the re- 
sulting jammed configurations are analyzed to determine 
whether they are collectively jammed and unique. 



III. ANALYSIS OF JAMMED STATES 

To verify if a given final configuration is collectively 
jammed we analyze the e igen value spectra of the dynam- 
ical (or rigidity) matrix |16| . Mi at jp, where the indices 
i and j refer to the particles and a, (3 = x, y represent 
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FIG. 5: The total number of distinct CJ states n^ * (circles) 
and the number of trials required to find them (squares) ver- 
sus system-size TV. The solid and dotted lines have slopes 
equal to 1.2 and 1.7, respectively. 

the Cartesian coordinates. For a system with Nf floaters 
and TV' = TV — Nf particles forming a connected network 
the indices i and j range from 1 to TV'. Thus, the dy- 
namical matrix has dN' rows and columns, where d = 2 
is the spatial dimension. By differentiating the interpar- 
ticle potential we find that the elements of the dynamical 
matrix with i j are given by |25| 

M-iaJfi — ^ija^ij /3 ) Cij^ija^ij (3 , (3) 

where = dV/drij and Cy = d 2 V/drfj, while those 
with i = j are given by 

M la ^ = -Y^M iaJ p. (4) 

j 

The dynamical matrix © and |Q| has dN' real eigen- 
values d of which are zero due to translational in- 
variance of the system. In a collectively jammed state 
no set of particle displacements is possible without creat- 
ing an overlapping configuration; therefore the dynamical 
matrix has exactly dN' — d nonzero eigenvalues. In our 
simulations we use the criterion > £ m j n for nonzero 
eigenvalues, where £ m i n = 10~ 6 is the noise threshold for 
our eigenvalue calculations. 

We note that our energy minimization algorithm for 
creating jammed states does occasionally yield a config- 
uration that is not collectively jammed. These states, 
however, are not considered in the current study. The 
number of trials that yield collectively jammed states out 
of the original n' t trials is denoted n t . The fraction of tri- 
als that give locally but not collectively jammed states 
(n' t — ni) /n' t decreases with increasing system size from 
« 5% at N = 6 to less than 1% for N > 12. 

We determine whether two collectively jammed states 
are distinct by comparing the sorted lists of the nonzero 
eigenvalues of their respective dynamical matrices. If the 
relative difference between two corresponding eigenvalues 
differs by more than £diff = 10 -3 , the configurations are 



TABLE I: Maximum number of trials performed n™ ax and 
fraction of CJ states obtained {n s /n!'° t )- mBX versus system size 
TV. 
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treated as distinct. By comparing the topology of the 
network of particle contacts in a representative sample of 
CJ states, we have found that this criterion is sufficient 
to reliably determine whether two states are distinct or 
identical. 

This procedure allows us to determine the number n s 
of distinct collectively jammed states at each fixed num- 
ber of independent trials n t . As expected, if two C J states 
have different packing fractions, they are distinct, with 
different contact networks and dynamical modes. This 
property holds with very high numerical precision — the 
packing-fraction difference of 10~ 13 already assures that 
the two states are distinct. 

However, it is not true that all collectively jammed 
states with the same packing fraction are identical. For 
example, the C J states shown in Fig. [2] have the same 
packing fraction, but they possess different contact net- 
works and eigenvalue spectra. This is a clear demonstra- 
tion that two collectively jammed configurations at the 
same packing fraction can have very different structural 
properties. 

We have also calculated the total number of contacts 
between particles TV C , i.e. the number of bonds that sat- 
isfy Tij < dij, in our slightly compressed jammed con- 
figurations. We find that the number of contacts in the 
collectively jammed states satisfies the relation ^3 Hil 

TV C > TV c min = 2 (dN' - d + 1) . (5) 

The minimum number of contacts required for mechan- 
ical stability of the system TV™ 111 can be calculated by 
equating the number of degrees of freedom to the number 
of constraints. Note that an extra constraint is required 
to prevent particle expansion. We have found that nearly 
all of the collectively jammed states have TV C = TV™ 111 ; 
fewer than 1% of these states have N c > TV" lm as shown 
in Fig. All configurations that are not collectively 
jammed have fewer contacts than TV™ 1 ". 

IV. RESULTS 

In the preceding two sections, we described our meth- 
ods for generating and counting distinct collectively 
jammed states. We will now present the results from 
these analyses. We will first discuss how the number 
of CJ states depends on parameters such as the num- 
ber of trials and system size. We then decompose the 
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FIG. 6: Probability distribution P(<f>) for obtaining a CJ 
state at <j> for AT = 10 (solid line) and N = 14 (dotted line) 
at (n s /n* ot ) max . The distributions at n 3 /n!; ot =0.05 (squares 
for N = 10 and triangles for N = 14) overlap those with 
larger n s /nl ot . 

probability density of obtaining a CJ state at a given 
packing fraction (Eq. (IJ) into the density p{<p) of CJ 
packing fractions and their frequency distribution 0(<p). 
We also consider under what conditions all of the possi- 
ble CJ states can be enumerated and determine whether 
strong conclusions can be made about the distributions 
of CJ states in large systems even though complete enu- 
meration is not possible. 

Our studies of the number of distinct CJ states n s 
versus the number of independent trials n t led to several 
surprising observations. First, we find that these systems 
possess a significant fraction of rare CJ states and thus 
an exponentially large number of trials are required to 
obtain nearly all states. Second, a master curve appears 
to describe n s (nt) for systems with N > 10, as shown 
in Fig. 21 (Each data point in this figure was obtained 
by averaging over at least 100 distinct permutations of 
the nt trials.) Our numerical results indicate that when 
n s is more than about 20% of the total number of dis- 
tinct CJ states n' ot , the curve n s (n t ) can be accurately 
approximated by 



= 1-A 



lo Sio -t^ 



(6) 



where A w 0.05. 

Our direct computations for small systems (N = 6, 8, 
and 10) and numerical fits to the master curve © for 
N = 12 and 14 indicate that both n\ ot and n* ot increase 
exponentially with system size as shown in Fig. [5] How- 
ever, both these quantities remain finite for any finite 
system. In particular, for the smallest system sizes, we 
increased the total number of trials by at least a factor of 
10 and did not find any new collectively jammed states. 
We also used several different algorithms for generating 
CJ states, e.g. compression and expansion of particles 
followed by relaxation using molecular dynamics with 
dissipative forces along fij and frictional forces perpen- 
dicular to fij , and these did not lead to any new C J states 
that were not already found using the protocol described 
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FIG. 7: Density of collectively jammed packing fractions 
p(4>) for (a) N = 10, (b) 12, and (c) 14 at n s /n* ot = 0.2 
(solid), 0.4 (dotted), 0.6 (dot-dashed), 0.8 (long-dashed), and 
1.0 (dashed). 

in Sec. [HI The maximum number of trials and fraction 
of CJ states obtained are provided in Tabled 

As indicated in Eq. IjTjI. the probability distribution 
P{4>) for obtaining a collectively jammed state at a par- 
ticular packing fraction tf> can be factorized into two com- 
posite functions: the density of CJ states p{<f>) and the 
frequency @(<p) with which these states occur. In our 
simulations, the distribution P{<p) is calculated from the 
relation 



P{<P)d<P 



n P { 



dip) - np(4>) 



(7) 



Here np(<p) is the total number of CJ states (counting 
all repetitions of the same state) with packing fractions 
below 4>. The density of CJ states is evaluated using an 
analogous relation 



p{4>)d4> 



n s (cp + dip) - n s (ip) 



(8) 



where n s {(p) is the number of distinct CJ states that have 
been detected in the packing-fraction range below <p. (In 
fact, we have used the number of distinct packing frac- 
tions to define p(ip) in place of the number of distinct CJ 
states n s ((f>). However, this does not affect our results 
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FIG. 8: Density of collectively jammed packing fractions 
p{4>) for TV = 8 (solid), 10 (dotted), 12 (dot-dashed), and 14 
(long-dashed) at (n a /n t s ot ) ln ^. 

because distinct states with the same 4 are rare in 2d 
bidisperse systems.) We note that both the probability 
density (0 and the density of CJ states {HJ are normal- 
ized to 1. The frequency distribution (3(4) = P(<j>)/p(<j>) 
is normalized accordingly. 

Below, we show how P(4), p(4>), and (3((j>) depend on 
the fraction of CJ states n s /n* ot and system size TV. To 
plot these distributions, we used 10 bins with the end- 
point of the final bin located at the largest CJ packing 
fraction <f> max for each TV. We recall that the distribu- 
tion of CJ packing fractions p(4) does not depend on the 
protocol used to generate the CJ states. The protocol 
dependence of the distribution P(4) IS captured by the 
frequency distribution (3(4)- 

The probability distribution P(4>) of CJ states is shown 
in Fig. for two small systems TV = 10 and 14. The 
results indicate that P(4>) depends very weakly on the 
fraction n s /nl ot of CJ states obtained — only 5% of the 
CJ states are required to capture accurately the shape of 
P(4>) for these systems. This result holds for all system 
sizes we studied, which implies that the distribution of 
CJ states can be measured reliably even in large systems 
|l6Lll9| . Note that the width and the location of the peak 
in P(4) do not change markedly over the narrow range 
of TV shown in Fig. El 

To see significant changes in P(4), the system size must 
be varied over a larger range. P(<j>) for TV = 18, 32, 64, 
and 256 are shown in Fig. ^ a t fixed number of trials 
n t = 10 4 . The width of the distribution narrows and 
the peak position shifts to larger <f> as the system size in- 
creases. In Ref. 16], we found that P{4>) for this 2d bidis- 
perse system becomes a (^-function located at 4o = 0.842 
in the infinite-system-size limit. What causes P(4) to 
narrow to a <5-function located at 4o when TV — > oo? Is 
the shape of the distribution P(4) determined primarily 
by the density of states p(<f>), or does the frequency dis- 
tribution (3(<t>) play a significant role in determining the 
width and location of the peak? We will shed light on 
these questions below. 

We first show results for p(4>) and (3(4>) as functions 
of the fraction n s /nl ot of distinct CJ states obtained. In 
Fig- HI p(4>) is shown for several small systems. In con- 
trast to the total distribution P(4>), the density of states 
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FIG. 9: Frequency distribution (3(4>) for (a) N = 10, (b) 12, 
and (c) 14 for n s /n!; ot = 0.2 (solid), 0.4 (dotted), 0.6 (dot- 
dashed), 0.8 (long-dashed), and 1.0 (dashed). 

p(4>) depends on n s /ri* ot significantly. For N = 10, a sys- 
tem for which we can calculate nearly all of the C J states, 
the curve p(4>) reaches its final height and width when 
n s /n* ot « 0.5. However, its shape still slowly evolves as 
n s /n* ot increases above 0.5; the low-(/> part of the curve 
increases while the high-</> side decreases. This implies 
that the rare CJ states are not uniformly distributed in 4>, 
but are more likely to occur at low packing fractions be- 
low the peak in p(4>)- Similar results for p(4>) as functions 
of n s /nl ot are found for N — 12 and 14. By comparing 
p(4>) at fixed n s /n t ° t , we also find that p(4>) narrows with 
increasing N. To further demonstrate that p(4>) narrows, 
the density of states is plotted in Fig.|S|for several system 
sizes at (n s /n^ ot ) max listed in Tabled 

The dependence of the frequency distribution (3(4) on 
the system size N and the fraction n s /n* ot of CJ states 
obtained is illustrated in Fig. EI Tbe results show that in 
contrast to the functions P(4>) and p(4), the distribution 
(3(4) achieves its maximal value at the highest packing 
fraction for which CJ states exist </> ma x- By comparing 
(3(4) for different system sizes at fixed n s /n' ot we find 
that 4max increases with increasing TV. 



The frequency 
strongly peaked at <p max as n s /n" s " u increases 
tion of (3(4) with n s /n* ot can be explained by noting that 



distribution (3(4) 
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FIG. 10: Frequency distribution /3n{4>) normalized by the 
peak value for (a) N = 10, (b) N = 12, and (c) N = 14 
at n s /nl ot = 0.2 (solid), 0.4 (dotted), 0.6 (dot-dashed), 0.8 
(long-dashed), and 1.0 (dashed). Least squares fits to expo- 
nential curves (thin solid lines) are also shown for the largest 
n 3 /nl ot lit each N. 

(3(<fr) = P{4>)/ f>{4>) an d that P{<j>) does not depend on 
n s /n* ot for n s /n* ot > 0.05 according to the results shown 
in Fig. El The density of states p(</>) and the frequency 
distribution f3(4>) must therefore behave in opposite ways 
to maintain constant P{4>)- As shown earlier in Fig. [7] 
the peak in p(<p) widens (for n s /n* ot < 0.5) and shifts to 
lower packing fractions as n s /n* ot increases. Thus, the 
distribution must decrease at low packing fractions 
and build up at large packing fractions with increasing 
n s /nf\ 

In Fig. we show the frequency distribution /3„ (<fi) = 
/3(0)//3max, which is normalized by the peak value /3 max - 
The results are plotted on a logarithmic scale. The fre- 
quency distribution varies strongly with </>; C J states with 
small packing fractions are rare and those with large 
packing fractions (4> « 0.83) occur frequently. We find 
that /?„(</>) is exponential over an expanding range of <j) 
as n s /nl ot increases. For N = 10, f3 n {4>) increases expo- 
nentially over nearly the entire range of <f> at n s /n* ot = 1. 
We see similar behavior for N = 12 and 14 in panels (b) 
and (c) of Fig. thus we expect (3 n {4>) to be exponen- 
tial as N s /Nl ot -> 1 for N > 10. We have calculated 
least-squares fits to 

n =A /3 exp(Bpil>) (9) 

for the largest rt s /n* ot at each system size. As pointed 
out above, the frequency distribution becomes steeper 
with increasing N; we find that Bp increases by a factor 
of 3.5 as N increases from 10 to 18 (not shown). Note 
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FIG. 11: Cumulative probability distribution Fh of CJ states 
in a narrow range of packing fractions d<f> for N=12. The 
index i denotes the position of the state in a list ordered 
by the frequency of occurrence and n s is the total number of 
states in the given interval. The solid lines correspond to bins 
centered on <j> = 0.73, 0.75, 0.77, and 0.79; the dashed lines 
labeled 1, 2, 3, and 4 correspond to bins centered on <j> = 0.65, 
0.81, 0.83, and 0.69, respectively. The width of each bin is 
A(/> = 0.02. The inset shows that the data are well-described 
by Eq. I10L where Af ~ 2.4 and a varies from 0.3 to 0.4. 

that reasonable estimates of Bp can be obtained even at 
fairly low values of n s /n* ot . 

We showed in Fig. ^3 that the frequency distribution 
is not uniform in <f>] in contrast, it increases exponen- 
tially with cj). Fig. 1111 shows another striking result; the 
frequency distribution is also highly nonuniform within a 
narrow range of 4>. In this figure, we plot the cumulative 
distribution of the probabilities of jammed states in 
an narrow interval dip versus the index i in a list of all 
distinct states in dej) ordered by the value of the probabil- 
ity of each state. The data for several different intervals 
appear to collapse onto a stretched exponential form, 

F h =exp[-A F (l-i/n s ) a ], (10) 

where n s is the number of distinct CJ states within d<p 
and the exponent a varies from 0.3 to 0.4. These results 
clearly demonstrate that CJ states can occur with very 
different frequencies even if they have similar packing 
fractions. 

From our studies of small systems, we find that both 
the density p(4>) of CJ packing fractions and the fre- 
quency distribution (3(4>) narrow and shift to larger pack- 
ing fractions as the system size increases. (See Figs. 
and E3) How do these changes in p(<f) and f3(4>) affect 
the total distribution P(<p) and can we determine which 
changes dominate in the large system limit? To shed 
some light on these questions, we consider the position 
of the peak in P((f>) with respect to the maximal packing 
fraction of CJ states </> ma x for several system sizes. In the 
absence of changes in p(4>) as a function of (f> — m axj the 
maximum of P(<f>) should shift toward <fi = </> ma x with 
increasing system size, because the frequency distribu- 
tion P becomes more sharply peaked at max according 
to the results in Fig. ^| However, as shown in Fig. ^] 
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FIG. 12: P ((/>') (thick lines) and ^ n ((/>') (thin lines) for 
N = 10 (solid), 12 (dotted), and N = 14 (dot-dashed) at 
(n s /n* ot ) max , where <£' = P(^) and ^(^) for AT = 10, 

12, and 14 have been shifted by = 0.013, 0.005, and re- 
spectively, so that /3max for the three system sizes coincide. 
j3 n (</>) for each N has also been amplified by a factor of ~ 20. 

we find the opposite behavior over the range of system 
sizes we considered: the peak of P(4>) shifts away from 
0max- This suggests that the density of states, not the 
frequency distribution, plays a larger role in determining 
the location of the peak in P(<f>) in these systems. 

Additional conclusions about the relative roles of the 
the density of states and the frequency distribution on 
the position and width of P(<fi) can be drawn from our 
observation that the frequency distribution [3 is an ex- 
ponential function of <fi (c.f. the discussion of results in 
Fig. 1 1 U|l and that P((j)) is Gaussian for sufficiently large 
systems (as shown in and illustrated in Fig. I13|) . If 
we assume that the exponential form of the frequency 
distribution @ remains valid in the large-system- limit, 
the density of states p(<j>) = P{<f))/ (3(<j>) is also Gaussian 
with the identical width <j(N). The location of the peak 
in P((f>) is 

^p(N) = ^* p (N) + B (N)a(N), (11) 

where <f>p(N)js the location of the peak in p(<fi). In previ- 
ous studies we found that the width of P(4>) scaled 
as a ~ N~ n , with ft w 0.55. We have also some indi- 
cation that Bpa 2 decreases with increasing system size: 
Bpa 2 = 0.017 at N = 12 compared to 0.012 at N = 18. 
However, we are not currently able to estimate Bp in the 
large N limit. 

If the system-size dependence of Bp is weaker than 
N 2n , the quantity Bpa 2 will tend to zero and the fre- 
quency distribution will not influence the location of the 
peak in P(<f>)- In this case P(<j>) becomes independent 
of the frequency distribution in the limit N — > oo for 
a class of protocols that are characterized by a similar 
frequency distribution [3 as our present protocol. Thus, 
as our preliminary results suggest, random close packing 
can be defined as the location of the peak in p((f>) when 
N — > oo, and this definition is completely independent of 
on the algorithm used to generate the CJ states. In the 
opposite case, where the system-size dependence of Bp 
is stronger than N 2n , the position of the peak in P(<j>) 
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FIG. 13: The distribution P(<f>) of CJ states for (a) N = 18, 
(b) 32, (c) 64, and (d) 256 are depicted using circles. The 
solid lines are least squares fits of the large-0 side of P{4>) to 
Gaussian distributions. 

results from a subtle interplay between the density of 
states and the frequency distribution. However, even in 
this case one can argue that the dependence of the posi- 
tion of the peak only weakly depends on the protocol: a 
shift of the peak position requires an exponential change 
in the frequency distribution (3{<p). 

V. CONCLUSIONS 

We have studied the possible collectively jammed con- 
figurations that occur in small 2d periodic systems com- 
posed of smooth purely repulsive bidisperse disks. The 
CJ states were created by successively compressing or 
expanding soft particles and minimizing the total energy 
at each step until the particles were just at contact. By 
studying small 2d systems, we were able to enumerate 
nearly all of the collectively jammed states at each system 
size and therefore decompose the probability distribution 
P{4>) for obtaining a CJ state at a particular packing frac- 
tion <p into the density p(<p) of CJ packing fractions and 
their frequency distribution @(<fi). The distribution P((j)) 
depends on the particular protocol used to generate the 
CJ configurations, while p(<fi) does not. This decompo- 
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sition allowed us to study how the protocol-independent 
p(<j>) and protocol-dependent &{4>) influence the shape of 

P{4>). 

These studies yielded many important and novel re- 
sults. First, the probability distribution P((f>) of CJ 
states is nearly independent of n s /n* ot , and thus it can 
be measured reliably even in large systems. This fin ding 
validates several previous measurements of P(4>) |!6lll9| . 
Second, the number of distinct CJ states grows exponen- 
tially with system size. In addition, a large fraction of 
these configurations are extremely rare and thus an ex- 
ponentially large number of trials are required to find all 
of the CJ states. Third, the frequency distribution (3((p) 
are nonuniform and increase exponentially with <fr. We 
also found that even over a narrow range of <f>, the fre- 
quency with which particular CJ states occur is strongly 
nonuniform and involves a large number of exponentially 
rare states. Finally, we have shown that P(<f>) becomes 
Gaussian in the large N limit. Since P(4>) — p(4>)(3(4>) 
and (3(<f)) is exponential, we expect that p{<j>) is also Gaus- 
sian and controls the width of P{4>) for large N. We also 
have preliminary results that suggest that the contribu- 
tion from /3(<j>) to the shift of the peak in P(4>) decreases 
with increasing N . We expect that p((f>) will determine 
the location of the peak in P{(j>) in the large N limit, 
and thus it is a robust protocol-independent definition of 
random close packing in this system. 



VI. FUTURE DIRECTIONS 

Several interesting questions have arisen from this work 
that will be addressed in our future studies. First, we 
have shown that the frequency with which CJ states oc- 
cur is highly nonuniform. It is important to ask whether 
the rare states can be neglected in analyses of static and 
dynamic properties of jammed and nearly jammed sys- 
tems. For example, we have shown that P(4>) is insensi- 
tive to the fraction n s /nl ot of CJ states obtained and thus 
P{4>) is not influenced by the rare CJ states. However, 
rare CJ states may be important in determining the dy- 
namical properties of jammed and glassy systems if these 
states are associated with 'passages' or 'channels' from 
one frequently occurring state to another. Moreover, an 
analysis of the density of states and the frequency dis- 
tribution (both as a function of <p and locally in (f>) may 
shed light on the phase space evolution of glassy systems 
during the aging process. 

A closely related question is what topological or geo- 
metrical features of configurational phase space give rise 
to the exponentially varying frequency distribution? Can 
one, for example, uniquely assign a specific volume in 
configurational space to each jammed state? A candi- 
date for such a quantity is the volume Q c of configuration 
space in which each point is connected by a continuous 
path without particle overlap to a particular CJ state. It 
is likely that those CJ states with large tt c will occur fre- 
quently for a typical compaction algorithm, while those 
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FIG. 14: P(<f>) (solid), p((j>) (dotted), and /3(<f>) (dot-dashed) 
for a 2d monodisperse system with N = 24. 

with small f2 c will be rare. 

Another important question is whether the results for 
p{4>), f3(4>), and P(<j>) found in 2d bidisperse systems also 
hold for other systems such as monodisperse systems in 
2d and 3d. Does p(4>) still control the behavior of P((f>) 
or does the frequency distribution play a more dominant 
role in determining P((f>)1 To begin to address these ques- 
tions, we have enumerated nearly all of the distinct col- 
lectively jammed states and calculated P(0), p(</>), and 
P(4>) in small 2d periodic cells containing N — 4 to 32 
equal-sized particles. 

In our preliminary studies, we have found several sig- 
nificant differences between 2d monodisperse and bidis- 
perse systems, which largely stem from the fact that par- 
tially ordered states occur frequently in the monodisperse 
systems. First, in 2d monodisperse systems there is an 
abundance of distinct CJ states that exist at the same 
packing fraction. For example, in a monodisperse sys- 
tems with N = 24, multiple distinct states occur at 19% 
of the CJ packing fractions compared to less than 1% in 
bidisperse systems with N = 14. Second, for the system 
sizes studied, quantitative features of the distributions of 
CJ states depend on whether N is even or odd. Third, 
P((p) can possess two strong peaks. For example, two 
peaks in P(4>) occur at <f>x rj 0.805 and 02 ~ 0.844 for 
N = 24 as shown in Fig. Moreover, the large-0 peak 
that corresponds to partially ordered configurations is a 
factor of three taller than the small-</> peak that corre- 
sponds to amorphous configurations. Finally, the max- 
imum in (3((f>) coincides with the large-0 peak in P{4>) 
and f}{4>) decays very rapidly as (p decreases. As shown 
in Fig. 1141 the rapid decay of f3(<p) significantly suppresses 
the contribution of the peak in p(cf>) to the total distri- 
bution P(4>). Thus, P(4>), which depends on the protocol 
used to generate the CJ states, may strongly influence 
the total distribution P(</>) even in moderately sized 2d 
monodisperse systems. 

Many open questions concerning monodisperse sys- 
tems in 2d will be answered in a forthcoming article 
|27). We will measure the shape of P{4>) as a function of 
system size and predict whether p(tfi) or /3(0) controls 
the width and location of the peak (or peaks) in the 
large N limit. The fact that f3(4>) strongly influences 
P(4>) at small and moderate system sizes explains why it 
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has been so difficult to determine random close packing 
in 2d monodisperse systems 01 — different protocols 
have yielded different values for <^ rcp [53, H3 ■ 
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